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Initial cracks present in an engineering structure 
may lead to premature failure and evaluation of any new design of 
major structure from this point of view is drawing ever 
increasing attention of mechainical engineering designers. A 
variety of techniques are in use to predict fracture behaviour of 
engineering structures starting from simple closed form solutions 
to hybrid finite element methods Eind now Finite Element 
Alternating Technique (abbreviated as FEAT). Versatility , 
efficiency and ability to solve a variety of " Real World 
problems are the major advantages of this technique. This 
technique is reported to be about ten times faster than Hybrid 
finite element method for the same accuracy . In the present work 
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this technique is applied to a three dimensional analysis of a 
finite thickness plate containing an embedded elliptical flaw. A 
three dimensional Finite Element Program is developed using 
twenty noded isoparametric elements. This together with Trefftz’s 
solution for a three dimensional problem to a loaded elliptical 
crack forms the basis of finite element alternating technique. 
Results obtained have been compared with available numerical 
solutions, and an excellent agreement is observed for the 
embedded flaw problem under consideration. 
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INTRODOCTIOH 

JLJL REYIEW ^ PREYIODS WORK 


Design of large engineering struct-ures such as 
submarine shells, large pressure vessels etc. involves a variety 
of design analyses to ensure reliable performance . An important 
part of this design analysis is fracture analysis to ensure 
performance in presence of initial manufacturing flaws of known 
sizes and flaws developed during service. This fracture analysis 
is primarily aimed at determining, acceptable loads for a given 
flaw size and shape and the also fatigue life as a function of 
flaw size and load. 

This is often achieved by solving the problem 
under static loads only. Problem of crack propagation path and 
crack growth rate do not gain considerable importance from design 
analysis point of view as onset of crack growth is considered as 
the end of useful life of the component under design. It is thus 
often sufficient to find out the loads at which crack starts 
moving. Such static analysis is becoming a routine check in 
defense equipment design, several automotive component design, 
pressure vessel designs in nuclear and chemical 
plants, and aerospace structural designs. 

It is therefore essential to have an efficient 
algorithm and software for carrying out fracture analysis for a 
variety of problems of practical importance. It is noted that 
practical situations often involve : 

1, State of stress which could be analysed using a three 
dimensional technique. 

2. surface or embedded crack configurations. Most of the 
applications mentioned in earlier pareigraphs are likely to have 
mainly surface or embedded flaws resulting from fatigue , 


material defects , welding etc. Application and use of through 
crack is considerably restricted from the viewpoint of real life 
situations. It is also realized (1,19,20) that actual crack 
geometries are often modeled by an elliptical or semi- elliptical 
surfaces . 

Thus keeping in view these practical stress states and 
crack configurations a three dimensional approach has become 
the necessity. With the criteria of practicality 2 ind efficiency 
as prime concern search for a suitable technique was made. Use of 
simple closed form solution is good for initial conceptual level 
analysis and for this purpose a virtual atlas of formulae (2) are 
available to solve a variety of problems. However, these often 
imply' over simplifications and closed form solutions are 
nonexistent for most real life structures; for ex 2 tmple, pipe 
junctions , forged drive train member etc. . Finite difference 
methods predates Finite element method aind were used by Chen et 
al. ( 3 ), Baker (4), and Shmuely et al. (5) for stress 
intensity factors in a finite thickness central crack tension 
plate. Wilkins (6,7) used this technique for complex fracture 
problem such as an internally or surface flawed cylinder 
subjected to sudden pressurization. These methods are not 
recommended due to enormous and prohibitively large computer time 
required though still used some times. 

There are three dimensional finite element methods 
using following approaches: 

1. fine mesh of constant strain solid brick elements 

2. Elements with -fr singularity , where r is the 
distance from the crack tip. 

3. Collapsed quarter point elements. 

Using densely packed three dimensional constant strain 
solid brick elements along the curved crack front by Miyamoto et 
al. (8) for proper modeling of the 1 / -fr singularity results in 
an inefficient use of computer time (9) . Computational 

efficiency is improved by incorporating 1 / -^r strain singularity 
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in displacement, formulation of elements. Such elements were used 
for Mode I problem by Newman (10) and multimode surface flaw 
problem in press’irized cylinders by Raju (11) .These elements 
required addition of special subroutines to a general three 
dimensional FEM code. However, use of collapsed quarter point 
isoparametric elements (12,13) did not require any such additions 
and commercial packages like SAP IV (14) were used successfully 
to solve the embedded crack problem. The maximum errors in 
numerical results, when using collapsed quarter point elements, 
when compared with the theoretical results (15) were 5% and 7% 
for a circular crack and an elliptical crack with an aspect ratio 
of 1.5 respectively. 

Displacement Hybrid Finite Element Method was 
used for three dimensional analysis of surface flaws in thick 
walled reactor pressure vessels and aerospace structural 
components by Atluri et al. (16,17) .This method is based on 
modified variational principle of the total potential energy such 
that exact asymptotic form of solution for singular stresses and 
strains in elements near the crack boundary is incorporated. The 
variational formulation leads to final algebraic equations which 
contain generalized global displacements and K stress intensity 
factors as unknown variable. Thus these equations upon solving, 
directly give Kj >KlI >Kxii » along with general 
displacements u,v,w at the crack front .This hybrid displacement 
FEM is computationally more efficient than those using quarter 
point collapsed elements. 

Finally, the Finite Element Alternating Technique 
(hence forth referred to as FEAT ) which has been employed in 
present work, originally developed by Kobayashi et al.(18) 
requires two analytical solutions as follows: 

(1) Solution to the problem of an infinite solid with am 
elliptical crack subjected to cubic normal loading on the crack 
surface. This limitation on the type of loading was due to the 
non-availability of a general solution for a general polynomial 
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type of normal loading on the crack face, and 

(2) Solution to a semi infinite body subjected to uniform normal 
and shear stress over a rectangular portion of the surface. 

In this work limitation of solution 1 was one of the 
reasons for causing about 10% error which was rather poor as 
compared with three dimensional Hybrid FEM (16) . In 1981 
Vi jay kumar and Atluri (19) derived a generalized solution 
procedure for an embedded elliptical crack subject to arbitrary 
crack-face tractions in an infinite solid. This procedure could 
be manually evaluated for a simple linear pressure variation of 
load on the crack face but required the use of a digital computer 
for a higher order pressure variation. And until 1983 when 
Nishioka et al. (20) evolved a generic procedure to numerically 
evaluate the required Elliptic Integrals in solution 1 , FEAT was 
considered as inaccurate though an efficient technique. In their 
approach, Nishioka replaced solution 1 and 2 of (18) by 

(1) Solution to the elliptical crack subjected to arbitrary 
normal and shear loading on the crack surface , in an infinite 
solid, and 

(2) Solution of uncracked body using a general ntnnerical 
solution technique such as finite element method. 

It was proved (19,20) that FEAT is not only a 
versatile technique in its application and capable of solving 
practical problems but also about 10 times faster than Hybrid 
Finite Element Method. Further this technique has been 
successfully used to solve multiple surface cracks problem in 
pressure vessels by O’Donoghue (1). 
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1^ OBJ ECTIVE AMP SgQPE OF PRESENT WORK 

The objective of present work is to determine stress 
intensity factors along elliptical crack boundary in a finite 
plate containing an embedded crack using three dimensional 
Finite Element Alternating Technique. Tensile loads are 
considered on the plate such that the problem is restricted to 
mode I analysis. 

In the present work loading on the crack face in 
solution 1 expressed as a fourth order polynomial. The work 
involves implementation of three dimensional Finite Element 
Method and solution 1 of FEAT in a single program . Twenty 
noded isoparametric elements have been used in FEM analysis and 
multiple linear regression is carried out in order to determine 
stresses at the surface of elements. 

Chapter II describes the FEAT technique as used in 
present work. Chapter III describes the program structure and 
related important features of program development for Finite 
Element Alternating Method .Also briefly explained are the 
variety of numerical techniques used in the program. 

Finally, in Chapter IV results have been included and 
discussed. Chapter V presents conclusions drawn and suggestions 
for further work. 
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£1A£1EB 11 

flHJTE £LSi^ ALTERSATIHG TECHNIQUE 

It is mentioned in chapter I that FEAT is oja. 
alternating use of FEM and a closed form solution for a loaded 
elliptical crack in an infinite body subjected to arbitrary crack 
face tractions. Therefore, before discussing FEAT in any further 
details the two solution procedures are described briefly. 

BU flHITE ELmgfil METHOD 

Since considerable amount of literature is available 
not only on the method (22) itself but also on programming the 
method (23), here only a brief treatment is included in order to 
highlight the salient features of the FEM as used in present 
work. The three dimensional Finite Element Method used has three 
displacement components as three degrees of freedom of each 
node. A twenty (20) noded isoparametric solid element is used in 
the present work. As FEM is being used for uncracked plate only 
there is no need of any singularity or special element in this 
method. This probably is the reason behind the high computational 
efficiency of FEAT. The displacement vector of an element is 
expressed in terms of the shape fxinctions and nodal 

displacement values as vinder 

20 

{ U } = 2 Ni Ui 

i=l 

20 

2 Ni Vi 

i=l 


{ V } = 


( 2 . 1 ) 
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20 

{ W } = 2 Ni Wi 

i=l 

where shape functions are expressed in terms of local 
coordinates S,T and U as under (22) 

for corner nodes with = +/- 1 , = +/- 1 , =+/-! 


Ni =(1+S0)(1+T0)(1+U0) (S0+T0+D0-2) / 8.0 (2.2) 

for midside nodes with = 0 , =+/-!, =+/-! 

Ni =(1-S^)(1+T0)(1+U0) / 4.0 (2.3) 

for midside nodes with Ti = 0 , Si =+/-!, Ui =+/-! 

Ni =(1-T^)(1+S0)(H-U0) / 4.0 (2.4) 

for midside nodes with Di = 0 , Si =+/-!, Ti =+/-! 

Ni+(1-U^)(l+S0)(l+T0) / 4.0 (2.5) 

where S 0 = S . Si 

T0 = T . Ti ' (2.6) 

U0 = U . Ui 


For calculation of stresses at the Gauss points nodal 
displacements and strain matrices are generally used ; however, 
for calculation of stresses at the surface of the elements as 
would be needed in FEAT this approach leads to highly inaccurate 
stresses and is inadequate. Therefore multiple linear regression 
analysis (24) has been used to calculate stresses on surface of 
the element . The stress component of interest <133 is expressed 
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as 

<^33 = a 0 + ai . XI + a 2 . X 2 + ^3 . X 3 (2.7) 

where xj^,X 2 ,X 3 are the global coordinates of the point at 
which stress is to be calculated. The Linear Regression Algorithm 
is used to find coefficients a 0 to 33 of multiple linear 
regression for the desired elements. 

It is important to mention here that Frontal solving 
technique for solving final simultcineous equations of FEM have 
been employed in the present work and the routine (23) used is 
suitable for only symmetric system of equations. Use of the 
algorithm by Nishioka, et al ( 20 ) for equation solving was 
developed by Mondkar 2 ind Powell (25) and is also a powerful out 
of core equation solver for FEM applications auad is based on 
Grout’s algorithm. Though no definite comparison between the two 
algorithms could be made, Frontal technique has been used in 
present work for following reasons 

1. Global Stiffness matrix does not have to be assembled 
in Frontal technique as in case of Nishioka’ s approach there by 
requiring less of core memory. 

2 . FRONT as used here has re-solution facility to 
greatly reduce solution time for subsequent iterations, which is 
claimed to be the main advantage of algorithm used by Nishioka 
( 20 ). 

Program developed interpolates coordinates of all 
midside nodes if they lie on a straight edge there by avoiding 
entry of coordinates of all nodes manually. It is to be noted 
that even though the input of coordinates of each node is a semi- 
automatic one , developing a new mesh does call for manual input 
of all connectivity and fixity data. 
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2 ^ SPLOTIOH EBQ CE P TOE £Qfi AH EItLI£IJi?AJi CBAgS 
IM AB J BFIBJ TE BODY SUBJECTED TQ 
ARBITRARY IIRACE FACE TRACTIONS . 


An generalized solution procedure for an elliptical 
crack in an infinite body subjected to arbitrary crack face 
tractions was developed by Vijaykumar and Atluri (19) . This 
solution is the second solution used in the finite element 
alternating technique. This solution is explained and relevant 
portions are derived in following paragraphs. The presentation of 
algebraic details of the solution procedure is however, kept to a 
minimum in the interest of clarity as well as reasons of space. 

This solution procedure is based on TREFFTZ 
FORMULATION (14) . Let Ux(x=l,2,3) and o^y (x,y=l,2,3) denote 
displacements and stresses, respectively in a homogeneous, 
isotropic linear elastic solid. The stress-displacement relations 
are, by Hooke’s law, in rectangular coordinates Xa(a=l > 2 , 3) 

2.V 

tTap = G (lict.p + Ua,(3 Sajj up^r ) (2.8) 

1 - 2.V 


where G and v are the shear modulus and Poisson’s ratio 
of the material , respectively, and 6^13 is dirac delta function. 
The Navier displacement equations of equilibrium in the absence 
of body forces are 

U(3,|3a + (1 - 2.V) Ua,|33 = ^ (2.9) 


In the foregoing the notation, ( ),a and ( ),ctp have 
been employed to denote first order and second order 
differentiation w.r.t. a and a,P respectively. 

Let R be the region of discontinuity in the plane X 3=0 
such that, after deformation, the material inside R breaks up 
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with free upper and lower surfaces, and remains continuous 
outside R. To deal with such a problem it is convenient to 
consider i+s complementary problem in which the surface of region 
of discontinuity is subjected to arbitrary tractions cscc- 

It is well known that (26) the solution for the 
problem can be expressed in terms of harmonic functions 8 and 
(a=l,2,3) in the form 

Ucc = + ^3 0>a ^2. 10) 

so that the equations (2.9) are satisfied if 

^a,a + ( 3 - 4v ) 8,3 = 0 (2.11) 

The stress components in terms of ® 

act|3 = G C ^a,|3 + ^f3,a + + ^1338,0 + 2x38, + 

2 v 

+ + 8,3 ) ] (2.12) 

(1 - 2 v) 

The problem is further simplified by expressing 8 and ^ 


in the form 

8 = ■7 .f = ' (2.13) 
<61 = (1 - 2 v)(fi ,3 + f 3 ,i) - (3 - 4 v)fi ,3 (2.14) 
^2 = (1 - 2v)(f2,3 + f3,2) - (3 - 4v)f2,3 (2.15) 
^3 = - (1 - 2 v)(fi,i + f 2 , 2 ) - 2(1 - v)f 3,3 (2.16) 


The stress components terms of f^ (cc— 1,2,3) are 

given by 
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ail = 2 G[f 3^11 2vf3^22 " 2fi^31 - 2vf2,32 + X 3 ( V . f ) ^ H] (2 . 17a) 

022 = 2G[f3^22 2vf3^11 - 2f2,32 “ 2vfi^31 + X 3 ( V ■ f ) , 22 ^ ( 2 . 17b) 

ai 2 = 2G[ ( l- 2 v)f 3 ^ 12 - (l-v)(fi^32 + f2,13> ^3 ( V . f ) , 12 ] 1 2 . 17c) 

033 = 2G[ - f 3,33 + X3(V.f).33] (2.17d) 

031 = 2G[-(1- v)fi ^33 +v(fi^ll + f2,21>+ X 3 (V.f),i 3 ] (2.17e) 

032 = 2G[-(1- v)f2,33 +v{fi,i 2 + f2,22)+ X3(V.f),23] (2.17f) 

where 

,1 ,2 and ,3 denote differentiation w.r.t. xi,X 2 

and X 3 respectively and 

V . f = fi,i + f 2,2 + ^3,3 (2.18) 

Now the problem is reduced to finding appropriate 
potential functions fa(ot = 1,2,3). Suppose that xi and X 2 are 
the cartesian coordinates in the plane of the elliptical crack 
Bind X 3 is normal to the crack-plane as shown in Fig. 2.1 such 
that : 

( xi/ai)^ + (X 2 /a 2 )^= 1» ai > a 2 (2.19) 

describes the border of the elliptical crack of aspect ratio 
(ai/a 2 ). The ellipsoidal coordinates $i{i=l,2,3) are more 
convenient in describing the above geometry and are defined as 
the roots of the cubic equation 


w(s) = 1 




. + s 


2 

X2 

2 . „ 
a2 + s 



= 0 


( 2 . 20 ) 


s 
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where Eq. (2.20) is written in an alternate form'- 

w(s)= P(s)/Q(s) ( 2 . 21 ) 

and 

P(s)= (s - $i)(s - $ 2 )(s - $ 3 ) (2.22a) 

Q(s)= s (s + ai^)(s + a 2 ^) (2.22b) 


where 


2 2 

w > Si 0 S $2 ^ " ^2 ^ $1 it - aj 

the elliptic boundary (Fig. 2.1) in the plane X 3 = 0 

corresponds to the curve $3 = $2 ~ crack surface itself 

namely the region inside ellipse in the plane X 3 = 0 is now given 
in a simple manner by the surface ^3 = 0 .The three potential 
functions fQ(a=l,2,3) of the Trefftz’s formulation are taken to 
be : 


fct = 2 2 Cc(,k,i Fki 

k 1 


(2.23) 


Fkl 


,k+l 


r 


I [w(s)3 

6x1^ 6x2^ J53 


k+ 1+1 


ds 


^Q(s) 


(2.24) 


In Eq. (2.23) are xmdetermined coefficients and the 

commas are introduced for convenience only . 
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FIG. 2-1 ELLIPTICAL CRACK IN AN INFINITE 
SOLID 
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Denoting by fa ,<3 partial derivatives of with 
respect to xp (0 = 1 , 2 , 3 ) , the first , second , and third partial 
derivatives of f^t are expressed as; 


fa,p =22 Ca,k,l ^kl,3 
k 1 


(2.25) 


fa, 136 = 2 2 Ccj.k, 1 5'kl,36 

k 1 


(2.26) 


fa,j36r =22 Ca,k,l Ffei.psr 

k 1 

By successive differentiation , it can be 
Eq . (2.24) that 


(2.27) 

shown f rom 


r' 


gk+i ^k+i+i 


r 


Fki = 1 ^ I 

J $3 6 xi^ 6 x 2 ^ *fQ(s) J $3 


ds 

klr-lltinl ^k+1+1 




w 


(2.28) 


4Q(s) 


wherein 


ki = k; li = 1 ; mi = 0 


(2.29) 


In Eq. (2.28) the additional notation that 6 ^*^ implies 
the jth partial derivatives w.r.t. x^ has been used. Similarly , 
the first order partial derivatives of Fki w.r.t. (a =1,2,3) 
cam be expressed by 


f 

Fkl,0 = |«3 


ds 


.klsll.ml „k+l+l 
6 6 6^^ w^ 


4Q(s) 


(2.30) 


where 


ki = k + 6 ip ; li = 1 + §20 ; mi = 63 ^ 


(2.31) 


In case of second and third partial derivatives, following 
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could be derived: 


f 

Fki,pr = Ua 

j 


Rkl.ll.ml k+1+1 
0, 6- 0- w 


ds 


4 Q(s) 



where k^ - k + 6^13 + 1^ = 1 + 623 + 62r; = 633+63^ 

6#v 


j* ds 

„ _ , ^klRll^ml k+1+1 . 

Fkl,pr6 - 1®3 Sj 63 w + 


^kipr 


+ G0 ( 2 . 33 ) 


J 


-fQCs) 6x5 


where kj = k0 + 615; Ij = Ig + 62; mi = in^ + 635. 


where k0 = k + dip + 6ir; I0 = 1 + 623 + ^ 2 T> ”“0 = ^30 S 3 r- 

0 

Numerical evaluation of terms G , Fi[i3r and their 
derivatives is vital to determination of stresses in the body 
due to crack loading . These numerical procedures are covered in 
detail in ( 20 ). Also it should be noted that for evaluation of 
second and third order derivatives of function Fj^i a generic 
integral of the following tsrpe needs to be frequently evaluated. 


f 


ds 


-kl r 11 .ml k+1+1 
S3 6, 6^ 63 w 


J 


WQ(s) 


( 2 . 34 ) 


This can be achieved by expanding in terms of x“^ and 

carrying out term by term differentiations. Then one gets 
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U3 

6“ 6“ 6”' 

= (k+l+l)l 

k+1+1 

2 

p q 

2 2 ( 

J 

"(0(3) 


p=0 

q=0 r=0 


(2p-2q) ! (2q-2r) ! (2r) ‘ 

^2p-2q-k ^2q- 

-2r-l 

^2r-m 


(k+l+l-p) ! (p-q) ! (q-r) ! (r) ! (2p-2q-ki) ! (2q-2r-li) ! (2r-mi) i 


■ ) (2.35) 

where ( . ) denotes multiplication and 

ds 

Jp-q,q-r,r(*3) = 1$3 (2.36) 

J {s+a*i)^ ^ (s+a®2)*^ ^ (s+a* 3 )^ •iQ{s) 

Thus one of the key algebraic steps in the successful 
application of FEAT is evaluation of generic elliptic integrals 
in equation (2.36). Nishioka et al (20) have developed a 
systematic procedure for the evaluation of this elliptic 
integral . 


EVALUATIQH QE 




UQ EIiLlfTIg IWTEqRMg 


In general the integrals ' in Eqs. 2.36 could be 
expressed in a closed form using incomplete elliptic integrals 
of the first and second typ^e and Jacobian elliptic functions, 
however, development of such expressions involves a lot of 
algebraic work and therefore there is a need for a systematic 
procedure for evaluation of these integrals as a function of 
parameter p, q, r. Thus rewriting equation (2.36) using Jacobian 
elliptic functions as 


'^p-q.q-r,r ^ 


2 


^2p^l 


tp>,q-r,r 


(2.37) 
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where 


"P,q-r,r - 




(2r-l)k’ = 


+ [2(-p+r-l) + 2(p-q-r+2) k* ]Lp^p_q^ j.-l 


+ k* ( 2p+2q~3 ) Lp ^ q— i* , r “2 ^ * 


(2.38) 


k» = (a»i - a*2) / a’l 


(2.39a) 


sn* uj^ = a*i / (a*i + $3) 


(2.39b) 


sn(u) ,nc(u) ,nd(u) are Jacobian elliptic functions. 

However, lowest order starting values of Lp^q-^^r-l 
and Lp^q-j;.,r-2 needed in order to evaluate Lp^q_j.,r- 

lowest order starting values are 


f 


ul 


^P,q,-1 ~ I0 (sn^^u) (nd^\i) (nc ^u) du 

J 


(2.40) 


ul 


r . ^ . 

^.q,-2 = I0 (sn^^u) (nd^Sa) (nc ^li) du 

J 

These integrals could be reduced ( 20 ) to 


(2.41) 


^p,q»-i ~ 


1 p 1 k>* 2 (i-r)p. 

2 2 


5-0 r=0 (p-j)!d! (i-r)!rs 


i2(q-j-r) 


( 2 . 42 ) 
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^p,q,-2 


1 


j,2p+4 


P 

2 

j=0 


2 j^,.2(2-r)^, 

2 

r=0 (P-J)!j! (2-r)!r! 


I2{q-j-r) 

(2.43) 


where 

|*Ul 

l2m = l0 u (2.44) 

J 

It is worth mentioning that in this series expansion of 
equation (2.42) denominator contains the term (l-Dl and not (2- 
F)! as given in (20). This has been verified by expanding 
integrals in Eq. (2.40 and 2.41) in terms of integrals 12® (27) 

of Eq. (2.44), 


Finally values of integral l2in various positive and 
negative values of m are to be determined. This can be done if we 
determine 


1. Inverse of Jacobian elliptic function sn in equation 
(2.39b) to determine uj . 

2. Incomplete elliptic integral of second kind E(uj^). 

Numerical procedure to invert Jacobian elliptic 
function and determine incomplete elliptic integrals of first and 
second kind are discussed in section 2.6. Remaining details for 
determining l2m could be found in (20,27). Thus in order to 
determine generic elliptic integral of Eq. (2.36) various arrays 
of integrals as outlined above need to be evaluated . 
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RELATIQH BETWEEN CRACK-FACE T R ACT I ONS 
POTENTIAL FDNCTIONS 


The t-ractions along the crack surface is expressed in 
the following polynomial form 


1 1 M 

O' 3q = 2 2 2 

i=0 j=0 m=0 


m 


_ ^2m-2n+i 2n+j 

^ ^ ,m-n,n * ^ 


n=0 


(a=l,2,3)(2.45) 


so that the values i,o specify the symmetries of the load 
with respect to the axes of th ellipse. The solution in terms of 
potential function is assumed in the form 


fa 


2 2 2 2 Ca,k-i,l F2k-21+i , 21+ j '' (a=1.2,3) (2.46) 

i=0 j=0 k=0 1=0 


The potential functions in Eq. (2.23) are to be \mderstood 
as that k and 1 replace (2k-21+i) and (21+j) respectively. It can 
be shown (20) that following linear algebraic equations result 
for a Mode I problem, between A coefficients of Eq. (2.45) and C 
coefficients of Eq. (2.46). 


1 

2 

i=0 


1 ^ 

2 {A3jiji-n,n^ ~ ^ 

i=0 i=0 


1 (_^jm+i+j+l ^ 1, 

2 { 2G 2 2 (-1)^ 

j=0 {2m-2n+i) ! (2n+j) ! k=m 1=0 

(2.47) 


[I^^’ Jli+1,L2,0(2>) 

m = 0, 1 , > M 


Cs'k-i.i 1: 


(2.48) 


n = 0, 1. m. 


Li = (m-n+k-l+i), L 2 = (n+l+J), 
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(2k+i+j+l)! (2Li+26ia)i {2L2+262a)5 

Ice (ot=1.2) (2.49) 

(k-m)! (L 2 + 62 a)i 

The relation between A and C coefficients in Eq. (2.47) 
can be summarized in a matrix form : 




{ A } = 

[ B 

3 { C } 



(2.50) 


Nxl 

NxN 

I Nxl 





where N 

is 

total 

number 

of coefficients A or 

C. For 

an 

incomplete 

polynomial 

, the 

number of terms depend i 

Qot only 

on 

the parameter 

M of Eq. 

(2.48) 

but also on parameters 

i and j 

of 


Eq. (2.45), 

STRESS JWTEMSXTY EACTOBS 

Once the coefficients C are determined by solving Eq. 
(2.50) for given loading on the crack surface, the stress 
intensity factors corresponding to these loads are computed from 
the following Eq. (19) for Mode I problem 


TT 

1 1 

M k 


jl/2 

A^'^^ I 2 

2 2 

(_2j2k+i+j 

aia2 

i=0 j=0. 

k=0 1=0 



1 COS0 


sin0 

(2k+i+j+l) ! 

( — - 

j2k-21+i^ 

j21+j 


aia2 ai 


^2 






C3,k-l,l (2.51) 


where 0 is elliptic angle and 


A = a“ 1 sin® 6 + a* 2 cos* 0 


(2.52) 
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For Mode II and Mode III equations for stress intensity 
factors could be found in (19,20). 




RICAL 


SY A Lg AT IOH QE IHCOMPLETE ELLIPTIC 
imCB M ig QE FIRST MD &1RB 


Since evaluation of incomplete elliptic integrals of 
first and second kind is necessary to elliptic integral 
evaluation of Eq. (2.36), the Arithmetic-Geometric Mean method 
and Landen’s Descending Transformation (28) are to be discussed . 
The parameter k or k’ , and crack size aj, a 2 are used to 
calculate parameters m and a by following (28) Eqs. ; 

m = k* ; a = sin ^(^m) ; (2.53) 


These parameters are now used to determine (a 0 ,b 0 ,C 0 ) 
to start the process of Arithmetic-Geometric Mean Method in which 
one proceeds to determine the number triples (a^.b^jC^), 
(a 2 .b 2 , C 2 ) j . . . ■ » according to following scheme of 
arithmetic and geometric means. 


30 = 1 

ai=(a0+b0)/2 

a2=(ai+bi)/2 


b 0 = cos a 

C 0 = sin a 

bi=(a0b0)^/^ 

ci=(a0-b0)/2 

b2=(aibi)^/2 

C2=(ai-bi)/2 


1/2 

aN=(aN-i+bN-l)/2 bN=(aN-ibN-l) CN=(aN-i-bN_i )/2 (2.55) 

This process is terminated at the Nth step when c^ 
becomes sufficiently small to the degree of accuracy required. In 

-7 

present work converging value of cj^ — 10 is taken. One can now 
determine incomplete elliptic integral of fist kind ie., F(ui) 
and of second kind ie., E(ui) using following relationships. 
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ui = F( ui ) = / (2"aN) (2.56) 

N 

E( ui ) = 2 Ci sin <t>^ + Ec. F( ui )/Kc (2.57) 

i=l 

where Kc = complete elliptic integral of first kind. 

Ec = complete elliptic integral of second kind and are 

given by 

Kc - X / ( 2 an ) (2.58) 

N 

Ec = Kc + Kc [ 1/2 2 2^ c*i 3 (2.59) 

i =0 

Angles in Eq. (2.56) and (2.57) are determined using 
following Landen’s descending Transformations. 

tan (^n+i - = (t>n/an) ^n» ^ 0 =sin ^(sn)(2.60) 

This procedure of inverting Jacobian elliptic function 
and determining desired elliptic integrals is implemented in the 
program. Particular caution is needed in implementing Eq. (2.60) 
since tan”^ ( tan (^h) ) could be considerably off the true value 
if proper correction is not applied 'to the computed value. 
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ZJL £I1?ITK ET.EMKWT ALTERHATIHG METHOD 

The finite element alternating technique as used in 
present work alternately uses following two solutions; 

1. Solution to the elliptical crack subjected to arbitrary normal 
and shear loading on the crack surface , in an infinite solid as 
discussed in sections 2.2 to 2.5. 

2. A general solution technique ie three dimensional finite 
element method as briefed in section 2.1. 

The basic concept of altematingly using these two 
solution procedures to give desired solution to the problem of a 
finite plate containing crack is presented in Fig. 2. 2. The 
solution of one method is used to find load for the other method 
for example FEM solution of an uncracked body is used to find 
stresses on the crack face, which are reversed and used as loads 
on the crack face of an infinite body. This is done in order to 
satisfy zero stress condition on the crack surface. Solution 1. 
is then used to determine stresses at the boundary of the finite 
body in consideration, these stresses are reversed and used as 
load for next FEM solution. The FEAT requires following steps as 
shown in Fig. 2.3. 

(1) Solve the uncracked body under the given external loads by 
using FEM. The tmcracked body has the seime geometry with the 
given problem except the crack. 

(2) Using finite element solution, stresses on the crack surface 
of original crack are computed. When the solution has converged 
stresses at the crack surface should be negligible, otherwise 
there would be some stresses still on the crack face. These 
stresses are termed residual stresses. Residual stresses in 
iteration two of this method are shown in Fig. 2.4. 
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FIG. 2-2 FEAT APPROACH; SUPERPOSITION OF VARIOUS CONFIGURATIONS 
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Fig. 2.3 Flow chart for finite element-alternating technique 
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(3) Compare the residual stresses in step ( 2 ) with a 
permissible stress magnitude. In the present study three percent 
of the maximum external applied stress is used for the 
permissible stress magnitude. 

(4) To satisfy the stress stress boundary condition, reverse the 
residual stresses. Then determine coefficients A of Eq. (2.45) by 
using following least squares method: 


f 

" 1 Sc 


. ip ,, 

- <J3a) 


dS (a 


1,2.3) 


(2.61) 


J 

where 03 ^ ' is the reversed residual stresses calculated by the 

(s) 

finite element solution, 03 ^ is defined by Eq. (2.45), Sc is the 
region of the crack, and Iq are the functionals to be minimized. 


Rewriting Eq. (2.45) in matrix form : 

(P T 

oza ^ i ^ { Aoc } (2.62) 

eind substituting Eq. (2.62) in Eq. (2.61), following relation 
between A coefficients «uid the residual stresses is obtained. 

{A}=[H3"^{Ra} (2.63) 

where 

^ T 

[ H ] = Isc J- > < L (2.64) 

J 

f 

{ Ra > = ISc t I* > ^3a dS (2.65) 

J 

(5) Determine the coefficients C of Eq. (2.46) in the potential 
function by solving Eq. (2.50). 

( 6 ) Calculate the stress intensity factors for the current 
iteration by substituting C coefficients in Eq. (2.51, 2.52). 
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(7) Calculate the residual stresses on the external surfaces of 
the body due to loads in step (4) using solution 1 ie Eq. (2.17a 
-2.17f). To satisfy the stress conditions of original problem , 
reverse the residual stresses and calculate equivalent nodal 
forces. This step is highly computer time intensive and takes 
considerable computational time. These nodal forces { Q } can be 
expressed in terms of coefficients C : 

{ Q }m = - C <5 3m { C } (2.66) 

and 

r 

t 3m = Ism C N [ n 3 [ p 3 dS (2.67) 

J 

where m denotes the nvimber of finite element, 

[ N 3 is the matrix of isoparametric element shape 
functions, 

[ n 3 is the matrix of normal direction cosines eind 
[ P 3 is the basis function matrix for stresses «ind can be 
derived from Eq. (2.17a - 2.17f). For the Mode I problem in hand 
this has been derived as*. 
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[ P 3 = 


M 

2 

k=0 

M 

2 

k=0 

M 

2 

k=0 

M 

2 

k=0 


k 

2 2G ( Fiji ^ ll+Fjci, 22+3^3 ^kl,31l) 

1=0 

k 

2 2G ( (1-2v)Fij 1^ 12+X3 Fki,3i2) 

1=0 

k 

2 2G X3 Fiji, 313 

1=0 ( 2 . 68 ) 
k 

2 2G (Fiji,22+2vFiji,ii+X3 Fki., 322 ) 

1=0 


M k 

2 2 2G X3 F]j1,323 

k=0 1=0 


M k 

2 2 2G (-Fij 1,33+X3 Fki,333) 

k=0 1=0 


In order to save computational time, the matrix [ G 3m 
for all elements on the boundary of the body is calculated prior 
to the starb of the iteration shovm in Fig. 2.2. 

(8) Consider the nodal forces in step (7) as external applied 
loads acting on the uncracked body. 

All steps in the iteration process are repeated vintil 
the residual stresses on the crack surface become negligible 
(step 3). Thus the solution for the given problem is the 
superposition of all the solutions obtained thus for and hence 
the stress intensity factors of all iterations are added to 
obtain desired solution. 



SSAflEB III 


FBOBtiEM FQRMOLATIOH MQ METHODOLOGY 


3.1 Introduction 

This chapter describes the details of Finite Element 
Alternating Technique as applied to the problem of a plate, 
containing 2 ua embedded crack, under unidirectional tension. A 
three dimensional technique needs to be used here as plate is of 
finite thickness. The program stmacture and techniques for 
minimizing CPO time have also been explained. 

3.2 Pro bl em description 

Stress intensity factors along the crack boundary of 
an embedded crack are determined using FEAT. The embedded crack 
is elliptic and is at the center of a finite sized plate (refer 
fig 3.1). 

The problem is analyzed using FEAT under the 
following assumptions. 

(i) Plate material is isotropic and linearly elastic. 

(ii) Small deformation theory is applicable for the purpose 
of calculating stresses using Finite Element Method. 

(iii) Plate is subjected to normal uniform stress loading on 

surfaces X 3 = ± L (fig 3.1), thus the stresses at the 

crack front are normal and hence all 

relations in FEAT are for Mode I. 
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3^ FINITE ELEMENT MESH 

The finite element discretisation of the uncracked body 
has considerable influence on convergence of the technique as 
well as accuracy of the solution. Fig. 3.2 shows the finite 
element mesh of the uncracked body. Full advantage is taken of 
reflective symmetry in three directions eind only one eighth of 
the plate has been modeled. Fig. 3.2 also shows that crack 
geometry is modeled by curved edged elements instead of straight 
edged elements as shown in Fig. 3.3, Though meinual input data 
preparation is relatively more involved with curved edged 
elements, following major adv 2 intage exists in using curved edged 
elements. In case of curved edged elements only six terms are 
needed in the polsoiomial of Eq. (2.45) for describing the load on 
the crack face, as against twenty one terms in case of straight 
edged elements (1) in order that satisfactory convergence is 
obtained. This implies that [ G matrices in Eq. (2.66) (which 
are calculated for each element on the external boundary of the 
body) , be 60 X 6 instead of 60 x 21 when straight edged 
elements are used. It can be said that for computation of these 
matrices if T CPU seconds are required by straight edged elements 
then only T/3 seconds or less would be needed by curved ones. 
Thus curved edged elements reduce the CPU time for computation 
of [Glm matrices to less than one third. 

3 ^ PROGRAM STBOCTSEE 

While developing the program attempts have been made to 
modularise entire calculations in subroutines and functions so 
that following objectives are fulfilled to the highest possible 
extent . 

(i) Unnecessary repetition of Calculations is avoided. This is 
very important to the efficiency of the program. 

(ii) Program build up is easy and each module could be separately 
tested and then assembled into the whole program. 
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FIG. 3-2 FEM MESH USING CURVED EDGED ELEMENTS 
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With this modularisation concept program was basically 
divided in two parts, one for FEM and the other for solution to 
the loaded crack problem referred to as FEAT. The FEM portion 
starts with subroutine FEM3D and all those called by FEM3D. This 
portion of the structure is given in Fig. 3. 4. It should be noted 
that coordinates of all the midside nodes lying on straight edges 
are interpolated by the program. 

Program structure of FEAT part of the program is shown 
in Fig. 3.5 and 3.6. Program structure clearly shows that 
suitable modules have been made that could easily be used and 
called by different subroutines. 

Computation of generic elliptic integrals of Eq. (2.36) 
is done for each of the 32 boundary elements, at each surface on 
boundary , at each of the 49 Gauss points and for each element of 
"the P matrix calculation ie. more than 57,600 times per element. 
Much higher number of factorial or Gamma function evaluation is 
done. It is important that these computations are done in least 
possible time. For this purpose Elliptic integral calculations 
are split in two subroutines. One subroutine precalculates 
Integral arrays of equations (2.38), (2.40) and (2.41) for 
different p, q, and r parameters for each Gauss point in 
consideration. Another subroutine now simply applies Eq. (2.37) 
using precalculated array of integrals to get the desired generic 
elliptic integral. 

It is worth noting that argument of factorial function 
does assume negative values in Eq (2.35) and for these, the 
corresponding function subroutine is not called at all instead 
following identity (28) for factorial of negative integer is 
taken . 


1 / n! = 0 ; n < 0. 


(3.1) 



MAIN 
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Fig. 3.6 Program structure of FEAT PART II 
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Computation of ellipsoidal coordinates 
ie. roots of polynomial of Eq (2.20) are repeatedly performed 
for each Gauss point for all boundary surfaces. This is done by 
the program using Graeffe’s root squaring method. 


3.5 Testing Elllgti. c Integral 

Evaluation of each module is significant and leads to 
considerable saving in program development. Testing of mainly 
modules relating to elliptic integral evaluation need to be 
briefed here. 

These subroutines have been tested for following crack 
related data and results have been compared with (20) in Fig. 3.7 
for parameters 1, m, and n given by 

ai =2, a2 = 1. as = 0, = -4, $2 = “1. *3 = 0-0625: 

l=p-q, m=q-r, n=r (3.2) 

It is clear from Fig. 3.7 that the values calculated by 
program match well with those of ref. (20). 



Elliptic integral, Jtmn 
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FIG. 3*7 ELLIPTIC INTEGRALS OBTAINED BY PROGRAM FEAT 
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BSSOL TS MU Discussion 


iJ. INTBODOCTIOn 

This chapter presents einalysis of the results obtained 
by using the FEAT progreun developed. Results are compared with 
the solution for the problem of finite thickness plate containing 
an embedded crack (20). 

The crack and plate geometry is as shown in Fig. 3.1. 
Plate thickness 2h to crack parameter 2a2 ratio of 0.75 is taken. 
The aspect ratio of crack ( a 2 / a^ ) is taken to be 0.4. Plate 
material is considered to have following properties 

Young’s modulus, E = 1000 Mpa 
Poisson’s ratio = v = 0.3 

FEM discretisation of the plate is as shown in Fig. 
3.2. The mesh consists of 80 elements and 516 nodes and 1545 
total degrees of freedom. There are four layers of 20 elements 
each. Plate size is also shown in Fig. 3.2. In the frontal 
technique Maximum FRONT WIDTH of -309 is encountered. 

In this mesh 32 elements are at the botmdary of the 
mesh (one eighth of the plate). For these elements [ G 3ni 
matrices are calculated at the start of the iteration process. 
Nximerical integration over each of the surfaces of these elements 
is done using 7x7 Gauss points as the terms of [ G 3 m involve 
polynomial of 13^^ degree. The loading on the crack face is 
approximated using six terms (1, x^ » 3^, x*x* , Xj , ) of a 
fourth order polynomial of Eq (2.45) . The parameter M, i, and 
j in Eq (2.45) for polynomial loading are taken to be M = 2; i = 
0, j = 0. In numerical integration of Eq (2.64) and (2.65) 5 x 5 
Gauss points were used for each element lying on crack face. 
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strfigp Intenslt-V Fact-nyR 


The met-hod has shown monotonic convergence and in four 
iterations the solution converged to the true answer. 
Convergence is checked using the normalized residual force on the 
crack face, which is given by 

f 

J Sc «^33 dS 


<70 7c ai a2 / 4.0 

The term in the denominator represents the force on the 
crack face, if it were subjected to the saune uniform distributed 
load as that acting on the plate end. In four iterations 
monotonic reduction of normalized residual forces on the crack 
face is shown in Fig. 4.1, Table 4.1. In last iteration this 
nondimensional , normalized force has a value of 0.028 which is 
satisfactory. This could be further lowered to 0.01 by changing 
the algorithm as would be discussed in next chapter. 

Fig. 4.2 , Table 4.2 shows convergence of stress 
intensity factors with successive iterations. Fig. 4.2 also shows 
variation of stress intensity factors Kj; along crack boundary 
w.r.t. angle 9 at the end of iterative process. It is importamt 
to note that stress intensity factors are higher at points nearer 
to the plate surfaces, ie. at angles 8 = 90“. This trend is 
similar to stress intensity factors at crack periphery of an 
elliptic crack in an infinite plate subjected to tension as shown 
in Table 4.2. It is only by normalizing the stress intensity 
factors of finite plate with those of an infinite plate, that a 
much clearer picture emerges. The variation of normalized stress 
intensity factors is shown in Fig. 4.3, Table 4.3. It is clear 
from the figure that the normalized stress intensity factors are 
higher at 0“ = 90“ , which leads to the conclusion that effect of 



Residual forces on crack face 


lA 


FIG. 4-1 



NORMALISED RESIDUAL FORCES 
ON CRACK FACE WITH SUCCESSIVE 
ITERATIONS 
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Table 4.1, Normalized Residual Force On Crack Face 


Iteration 

Number 


Normalized Residual 
Force on crack 
face 


0.997 


0.640 


0.630 


4 


0.028 


Stress intensity factor, Ki,MPa v/Tn 



FlG.4-2 CONVERGENCE OF STRESS INTENSITY FACTORS 
WITH SUCCESSIVE ITERATIONS 
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Table 4.2, Stress in.tensi't 7 factors 


Elliptical Stress intensity factors Kj, MPa 


angle 

0“ Iter 1 

Iter 2 

Iter 3 

Iter 4 

Infinite Plate 

0 

30.78 

48.90 

66.51 

67.29 

61.62 

10 

31.93 

50.82 

69.16 

69.98 

64.52 

20 

34.69 

55.47 

75.64 

76.54 

69.46 

30 

37.96 

61.10 

83.58 

84.58 

76.00 

40 

41.07' 

66 . 66 

91.51 

92.62 

82.21 

50 

43.75 

71.62 

98.68 

99.89 

87.56 

60 

45.88 

75.70 

104.68 

'105.96 

91.83 

70 

47.42 

78.74 

109 . 14 

110.50 

94.93 

80 

48.35 

80.60 

111.92 

113.32 

96.78 

90 

48.66 

81.23 

112.85 

114.27 

97.40 


Normalised stress intensity factors, 
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0 10 20 30 40 50 60 70 80 90 

Elliptical angle 0 (deg.) 

FIG. 4-3 NORMALISED STRESS INTENSITY FACTORS, A 
COMPARISON 
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Table 4.3: Normalized Stress Intensity Factors 


Elliptic Angle 

9" 

Normalized 

Present 

• 

Stress 

Work 

Intensity Factors 

Nishioka (20) 

0 

1.092 


1.10 

10 

1.095 


1.10 

20 

1.102 


1.105 

30 

1.113 


1.12 

40 

1.127 


1.14 

50 

1.141 


1.16 

60 

1.154 


1.18 

70 

1.164 


1.20 

80 

1.171 


1.208 

90 

1.173 


1.21 


finite plate thickness of the plate is more at the points, nearer 
to the plate surfaces. Fig. 4.3 and Table 4.3, also compare 
normalized stress intensity factors obtained in present work with 
those of Nishioka et al.(20). Since the maximum error of 3.1 % 
occurs at 0 = 90® , it cein be said that results compare well with 
those of Nishioka. The error could however be due to use of only 
linear regression analysis in step (2) of FEAT (Fig. 2.3). 
Nishioka may have instead used Multiple polynomial regression for 
this purpose. 

The program was developed and finally run on ND-500 
computer. In the present work with the 80 elements, 1545 degrees 
of freedom mesh a total of 1530 CPD seconds were required for 
final convergence. Almost one third of this ie. 500 seconds are 
needed for computation of [ G ] m matrices for all elements on 
the boiindary of the mesh. 





5 JL OfiPcl H s jlQng 

Need of an efficient, algorithm to solve problems of 
practical importance is evident from the discussion in Chapter I . 
In the present work stress intensity factors at the crack front 
of am elliptical embedded flaw, in a finite thickness plate was 
solved using Finite Element Alternating Technique. Based on this 
work following conclusions could be drawn. 

1. A medium mesh is sufficient for convergence to the true 
answer and a fine mesh is not required at all. 

2. The technique is about ten times faster than Hybrid Finite 
Element Method. 

SigflPS ffgxk 

Following further work is suggested in order that more 
complicated component eind crack geometries could be handled. 

1. Use of multiple polsmomial regression in step (2) of FEAT. 

2. Programming equations relating to surface crack and 
multimode applications. 

3. Coupling of the Trefftz’s solution code developed, to a 
general FEM code e.g. SAP IV, or SESAM. 
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